/*
* INSANE - Interactive Structural Analysis Environment
*
* Copyright (C) 2003-2005
* Universidade Federal de Minas Gerais
* Escola de Engenharia
* Departamento de Engenharia de Estruturas
* 
* Author's email :    insane@dees.ufmg.br
* Author's website :  http://www.dees.ufmg.br/insane
* 
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* as published by the Free Software Foundation; either version 2
* of the License, or any later version.
* 
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
* GNU General Public License for more details.
* 
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA  02111-1307, USA.
*/
package br.ufmg.dees.insane.model.femModel.problemDriver;

import java.util.ListIterator;

import br.ufmg.dees.insane.analysisModel.AnalysisModel;
import br.ufmg.dees.insane.materialMedia.degeneration.Degeneration;
import br.ufmg.dees.insane.model.femModel.element.Element;
import br.ufmg.dees.insane.model.femModel.element.ParametricElement;
import br.ufmg.dees.insane.shape.Shape;
import br.ufmg.dees.insane.shape.ShapeC0;
import br.ufmg.dees.insane.util.IMatrix;
import br.ufmg.dees.insane.util.IVector;

/** 
 * A class representing the technique of selective integration
 * used to prevent the "locking" in the analysis of thin plates
 * of Reissner-Mindlin.
 *
 * @author Samir
 * @since 20/07/2006
 */
public class RMThinPlateSubstituteShearStrain extends ParametricPhysicallyNonLinearSolidMech {

	private static final long serialVersionUID = 1L;

	/* (non-Javadoc)
	 * @see br.ufmg.dees.insane.model.femModel.problemDriver.ProblemDriver#getC(br.ufmg.dees.insane.model.femModel.element.Element)
	 */
	public IMatrix getC(Element e) {
	    int numberOfDegreesOfFreedom = e.getNumberOfDegreesOfFreedom();
		IMatrix stiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    stiffnessMatrix.setZero();
		IMatrix finalStiffnessMatrix = new IMatrix(numberOfDegreesOfFreedom,numberOfDegreesOfFreedom);
	    finalStiffnessMatrix.setZero();
		ListIterator gns = ((ParametricElement)e).getDegenerations().listIterator();
	    AnalysisModel analysisModel = e.getAnalysisModel();
	    ShapeC0 shape = (ShapeC0)e.getShape();    
	    IMatrix cartesianNodalCoordsMatrix = e.getCartesianNodalCoordsMatrix();
	    IMatrix samplingPoints = shape.samplingPoints();
		IMatrix C = this.C(shape,samplingPoints,cartesianNodalCoordsMatrix);
		IMatrix Bc = this.Bc(analysisModel,shape,samplingPoints,cartesianNodalCoordsMatrix);
		while(gns.hasNext()) {
		    Degeneration degeneration = (Degeneration) gns.next();
	        double[] naturalCoords = degeneration.getRepresentation().getNaturalcoords().toDouble();
	        IVector shapeFunction = shape.getShapeFunction(naturalCoords, cartesianNodalCoordsMatrix);
	        IMatrix derivedShapeFunction = shape.getDerivedShapeFunction(naturalCoords, cartesianNodalCoordsMatrix);
	        IMatrix secondDerivedShapeFunction = shape.getSecondDerivedShapeFunction(naturalCoords, cartesianNodalCoordsMatrix);
			IMatrix matrixB = analysisModel.getInternalVariableOperator(derivedShapeFunction, secondDerivedShapeFunction, shapeFunction, cartesianNodalCoordsMatrix);
			IMatrix J = new IMatrix(2,2);			J.mul(derivedShapeFunction,cartesianNodalCoordsMatrix);
			IMatrix APT = shape.getAPT(naturalCoords,cartesianNodalCoordsMatrix);
			matrixB = this.modifyInternalVariableOperator(matrixB,this.getSubstituteMatrix(J,APT,C,Bc));
			IMatrix matrixC = degeneration.mountC();
			IMatrix k1= new IMatrix(numberOfDegreesOfFreedom, matrixC.getNumCol());
			k1.mulTransposeLeft(matrixB, matrixC);
			stiffnessMatrix.mul(k1,matrixB);
			double weight = degeneration.getRepresentation().getWeight();
			double jac = analysisModel.getJacobianTransformation(derivedShapeFunction,cartesianNodalCoordsMatrix);
			double coeff = e.getElmCoefficient();
			stiffnessMatrix.setScale(weight*jac*coeff);
			IMatrix matrixT = analysisModel.getTransformationMatrix(derivedShapeFunction,cartesianNodalCoordsMatrix);
			stiffnessMatrix.mul(matrixT);
			matrixT.transpose();
			matrixT.mul(stiffnessMatrix);
			finalStiffnessMatrix.add(matrixT);
		}
		return(finalStiffnessMatrix);
	}
	
	private IMatrix C(ShapeC0 shp, IMatrix sp, IMatrix cn){
		IMatrix c = new IMatrix(2*sp.getNumRow(),2*sp.getNumRow());
		c.setZero();
		IMatrix jacobian = new IMatrix(2,2);
		double [] samplingPoint = new double[2];
		int k = 0;
		for (int i = 0; i < sp.getNumRow(); i++) {
			jacobian.setZero();
			for (int j = 0; j < sp.getNumCol(); j++) {
				samplingPoint[j] = sp.getElement(i,j);
			}
			jacobian.mul(shp.getDerivedShapeFunction(samplingPoint,cn),cn);
			for (int J = 0; J < jacobian.getNumRow(); J++) {
				for (int K = 0; K < jacobian.getNumCol(); K++) {
					c.setElement(J+k,K+k,jacobian.getElement(J,K));
				}
			}
			k = k + 2;
		}
		
		return c;
	}
	
	private IMatrix Bc(AnalysisModel am, Shape shp, IMatrix sp, IMatrix cn){
		IMatrix bc = new IMatrix(2*sp.getNumRow(),cn.getNumRow()*3);
		bc.setZero();
		double [] samplingPoint = new double[2];
		int k = 0;
		for (int i = 0; i < sp.getNumRow(); i++) {
			for (int j = 0; j < sp.getNumCol(); j++) {
				samplingPoint[j] = sp.getElement(i,j);
			}
			IMatrix Bc = this.getInternalVariableOperatorShear(am.getInternalVariableOperator(shp.getDerivedShapeFunction(samplingPoint,cn),shp.getSecondDerivedShapeFunction(samplingPoint,cn),shp.getShapeFunction(samplingPoint,cn),cn));
			for (int J = 0; J < Bc.getNumRow(); J++) {
				for (int K = 0; K < Bc.getNumCol(); K++) {
					bc.setElement(J+k,K,Bc.getElement(J,K));
				}
			}
			k = k + 2;
		}
		
		return bc;
	}
    
	/**
	 * Mounts the internal variable operator shear matrix (BS).
	 * @param B matrix complete 5 x 5.
	 * @return bShear matrix, part of B matrix complete relative the shear.
	 */
    private IMatrix getInternalVariableOperatorShear(IMatrix B){
    	IMatrix bShear = new IMatrix(2,B.getNumCol());
    	bShear.setZero();
    	
		for(int i = 0; i < B.getNumCol(); i++){
			bShear.setElement(0,i,B.getElement(3,i));
			bShear.setElement(1,i,B.getElement(4,i));
		}
		return bShear;
    }
    
    private IMatrix getSubstituteMatrix(IMatrix J, IMatrix APT, IMatrix C, IMatrix Bc){
    	IMatrix bc = new IMatrix(2,Bc.getNumCol());    	bc.setZero();
    	IMatrix aux1 = new IMatrix(2,APT.getNumCol());    	aux1.setZero();
    	IMatrix aux2 = new IMatrix(2,C.getNumCol());    	aux2.setZero();
    	IMatrix invJ = new IMatrix(J.getNumRow(),J.getNumCol());	invJ.setZero();
    	
    	invJ = J.pseudoInverse();
    	aux1.mul(invJ,APT);
    	aux2.mul(aux1,C);
    	bc.mul(aux2,Bc);
    	
    	return bc;
    }
    
    private IMatrix modifyInternalVariableOperator(IMatrix B, IMatrix Bs){
    	IMatrix b = new IMatrix(B.getNumRow(),B.getNumCol());
    	b.setZero();
		for(int i = 0; i < B.getNumCol(); i++){
			b.setElement(0,i,B.getElement(0,i));
			b.setElement(1,i,B.getElement(1,i));
			b.setElement(2,i,B.getElement(2,i));
			b.setElement(3,i,Bs.getElement(0,i));
			b.setElement(4,i,Bs.getElement(1,i));
		}
    
    	return b;
    }

}